function pass = circplate
% Use a linear chebop to compute the deflection of a circular plate.

tol = chebfunpref('eps');

d = [0 1];
r = chebfun('r',d);
f = exp(-50*r);                    % loading function

B = chebop(d);                     % r^3*(biharm)
B.op = @(r,u) r.^3.*diff(u,4) + 2*r.^2.*diff(u,3) - r.*diff(u,2) + diff(u);
B.lbc = @(u) [diff(u), diff(u,3)]; % symmetry at origin
B.rbc = @(u) [u, diff(u,1)];       % clamped at boundary

u1 = B\(r.^3.*f);
u2 = chebfun(solution,d);

pass = norm(u1-u2) < 2e-9*(tol/eps);

function vals = solution
% Previously computed solution!
vals = 1.0e-04 * [0.492425828933118
   0.492420508896867
   0.492341044411885
   0.492000637352868
   0.491105572455726
   0.489291166029874
   0.486172420678670
   0.481393480475663
   0.474661932485259
   0.465764953929333
   0.454572917922542
   0.441037402172374
   0.425187532682987
   0.407125584553868
   0.387021455334826
   0.365105575259274
   0.341660168540615
   0.317009078624271
   0.291506527986001
   0.265525242144417
   0.239444376330276
   0.213637667843278
   0.188462208017432
   0.164248188715530
   0.141289931149257
   0.119838450980205
   0.100095754592482
   0.082210998782222
   0.066278581633689
   0.052338167847845
   0.040376589050578
   0.030331500365460
   0.022096620406903
   0.015528334291951
   0.010453399541128
   0.006677463865837
   0.003994082586928
   0.002193912281901
   0.001073756411306
   0.000445148028440
   0.000142173840725
   0.000028272198672
   0.000001774114372
  -0.000000000000000];
